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^ ; Abstract 

^ ■ Numerical solution of equations governing time domain simulations in com- 

•^ . putational electromagnetics, is usually based on grid methods in space and 

on explicit schemes for the time evolution. A predefined grid in the problem 
domain and a stability step size restriction must be accepted. Evidence is 
given that efforts need for overcoming these heavy constraints. Recently, the 
authors developed a meshless method to avoid the connective laws among 
the points scattered in the problem domain. Despite the good spatial proper- 
^ ! ties, the numerical explicit integration used in the original formulation of the 

method provides, also in a meshless context, spatial and time discretization 
strictly interleaved and mutually conditioned. Afterwards, in this paper the 
t^ ' stability condition is firstly addressed in a general way by allowing the time 

'nI" . step increment get away from the minimum points spacing. Meanwhile, a 

^^ '. formulation of the alternating direction implicit scheme for the evolution in 

time is combined with the meshless solver. The formulation preserves the 
leapfrog marching on in time of the explicit integration scheme. The new 
k>i ' method, not constrained by a gridding in space and unconditionally stable 

?H ■ in time, is numerical assessed by different numerical simulations. Perfect 

- - -' matching layer technique is used in simulating open spatial problems; oth- 

erwise, a consistency restoring approach is introduced in treating truncation 
at finite boundary and irregular points distribution. Three case studies are 
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investigated by achieving a satisfactory agreement comparing both numerical 
and analytical results. 

Keywords: ADI leapfrog method, meshless methods, smoothed particles 
electromagnetics. 



1. Introduction 

Meshless methods have recently emer ged as numerical techniques for elec- 
tromagnetic modelling [H, S, 0, H, 0, S, llil Q [isl M, 0, H- These 
methods do not require a predefined mesh, and use points scattered in the 
problem domain avoiding the need of information on the position among 
them. One of the most popular meshless method, the smoothed particle 
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hydrodynamics (SPH) |17|, ll8|, ll9|, l20|, |2l|, [221 123|, |2j has been recently refor- 



mulated by the authors for computational electromagnetics (GEM) problems. 
Maxwell's equations, which relate the electric and magnetic fields by means 
of a time dependent system of partial differential equations(PDEs) are con- 
sidered [l|, 121 . Due to the coupling nature of the electric and magnetic field 
components, the points are in such a way that the electric points should be 
surrounded by those of magnetic field and vice-versa. The method named as 
SPEM is applied in approximating the space field variables of the time do- 
main Maxwells curl equations by using a kernel representation working with 
a cluster of scattered nodes. The differential operators applied to the field 
variables are transmitted to a kernel function which characterizes the spatial 
discretization. Despite the good spatial properties of SPEM, the numerical 
integration based on explicit method for time evolution, introduces a severe 
constraint. In fact, as for the conventional grid based methods, SPEM has 
to satisfy the Courant-Friedrich-Levy stability condition [g, |lO|, |28| . As well 
known, this condition becomes highly restrictive when the geometry resolu- 
tion forces the use of much finer space discretization than that is needed to 
solve for waves away from material irregularities or interfaces. This problem 
is firstly investigated in the paper in SPEM framework, in order to deter- 
mine the largest allowed time increment. Moreover, in technical literature, a 
great interest is addressed in finding numerical schemes which allow cornpu- 
tational advantages together with unconditional stability |8|, |25|, |29|, |30|, |31 . 
Otherwise, an alternate implicit direction method is proposed for the time 
evolution of the electromagnetic field in the SPEM framework by preserv- 
ing the leapfrog marching on in time of the explicit integration scheme. By 



considering a 2-stage splitting approach, the leapfrog implicit formulation is 
carried out by suitably handling the formulas related to two adjacent time 
steps. By working with an implicit scheme, a sparse linear system is gener- 
ated and it has to be solved at each time step. The system matrix involves 
the kernel derivatives, and it can be assembled only once at the beginning 
of time integration process. Details on the matrices used in the problem are 
investigated. The resulting algorithm is called as LAF-SPEM method. Ac- 
curacy degradation, due to a lack of particle consistency, is also taken into 
account: a numerical corrective strategy, which allows to restore the con- 
sistency, without any modification of the smoothing kernel function and its 
derivatives, is employed [l| . Test problems are simulated to validate the pro- 
posed methodology. Simulations running with a time step larger than that 
satisfying the CFL condition has been found to remain stable, so improving 
the conventional SPEM algorithm. The paper is organized as follows. In 
section 2 basic ideas on smoothed particle hydrodynamics method are re- 
called in order to provide the original SPEM numerical scheme in solving 
time-domain Maxwells curl equations. Investigations about the stability cri- 
terion are also included. In section 3 the LAF-SPEM method is described 
and simulation results are reported and discussed in section 4 also by com- 
paring with standard FDTD computations and existing exact solution. A 
brief conclusion complete the paper. 

2. The original smoothed particle electromagnetics method 

The SPH scheme has been adopted by the authors to compute the time- 
domain Maxwells curl equations by using particle approximation, so obtain- 
ing SPEM method |l|, |2|. For the reader convenience, a brief summary of 
SPEM fundamentals is reported in the following. Let consider the standard 
Maxwell ' s curl equations in time domain, in a linear medium in the compact 
curl notation: 

V X H = J + eAE (1) 

V X E = -/iAH 

where E and H are the electric and magnetic vectors fields, J = aE is the 
conduction current density vector, t as the time, Dt the derivative operator, 
e, /i, a the medium permittivity, permeability, and conductivity, respectively. 



which are supposed to be constant. The SPH method is based on an integral 
representation of the unknown function u(x): 



u{x) = / u{y)6{x - y)dy (2) 

Jn 

where O C K'' is the problem domain and 5{-) is the Dirac delta function. 
The delta function is approximated by a smoothing kernel function W{-) 
depending on the spatial coordinates and on the smoothing length parameter 
h. The smoothing kernel function must verify the following properties: 

a. lim/,^0 W{x -y,h) = 6{x - y); 

b. W{x — y,h) > on a sub-domain of fi and W{x — y,h) = outside; 
c- InW{x-y,h)dy = 1; 

d. W{x — y, h) function value for a particle is monotonically decreasing 
with increasing distance away from the particle. 

The discrete spatial formulation of equations ([1]) has been obtained by us- 
ing SPH particle approximation: the integral representation is approximated 
over a set of particles placed in the problem domain. The particles are the 
points in which the fields components are computed at each time step, by 
using the information belonging to the neighbouring ones. The smoot hing 



length h defines the size of the support domain of each particle [13|, Il8|, [19 
. Particular attention has to be paid on the choice of h. A too small h can 
lead to an inaccurate solution because not enough particles may be placed 
inside the influence area of a fixed particle Xj ; a too large smoothing length, 
can smooth out local properties and the particle approximation suffers too. 
From now on, for each referred fixed particle Xi , an own smoothing space 
length hi will be considered. The SPH particle approximation of ([2]) can be 
expressed as a linear combination of translates of a kernel basis function: 



NP 

'<^i^i) = J2 u{xj)W{xi - Xj, hi)AVj (3) 

where NP is the total number of particles in the influence domain of the fixed 
one Xi and AVj measures the surrounding media at the position of particle Xj. 
NP is strictly related to the choice of the smoothing length hi. In computing 
the two interlaved Maxwell's curl equations, particles for the electric field 



(i?-nodes) and magnetic field (if-nodes) must be considered in the spatial 
domain. Due to the electric and magnetic coupling nature of Maxwell' s 
equations, spatial staggering of E-nodes and i7-nodes is required. The E- 
nodes and the i7-nodes are arranged so that each E-node is surrounded 
by i7-nodes and each if -node is surrounded by E-nodes. The electric field 
update depends on the magnetic particles in the neighbourhood of the fixed 
electric particle and viceversa. In solving ([T]), fields spatial derivatives have 
to be approximated by means of SPH. The differential operators on the fields 
components are transmitted to the smoothing kernel function: so the kernel 
approximation allows differential operations to be determined from the values 
of the function and the derivatives of the kernel, rather than the derivatives 
of the function itself. However, this is really true only when the influence 
domain of a flxed particle is placed entirely inside the problem domain; when 
the influence domain overlaps with the geometry boundary the smoothing 
kernel W is truncated, and a resulting non-zero surface integral has to be 
opportunely treated in order to avoid the corruption of the numerical solution 



18l |. In the following, for the sake of simplicity, and in agreement with the 
numerical results reported in section 4, a lossless medium is considered. Thus, 
the two interleaved Maxwell's curl equations (j2]) result in SPEM formulation 



(4) 



as follows: 




DtH.irf) 


NPh 

= -I E [E.{rf)DyW{rf - rf , hf)- 




Ey{rf)D,W{rf-rf,h^)]AV, 


DtHy{Tf) 


NPh 

= -^ E [E,{rf)D,W{r^-rf,hf)- 



E,{rf)D,W{rf - rf, h^)]AVj 



NPh 

D.H^irf) = -^ E [Eyirf)D^Wirf - rf, hf) 
E,{rf)DyW{r^ - rf, hf)]AVj 



NPe 

DtE^rf) = i E [H.{rf)DyW{rf - rf, hf)- 
Hy{rf)D^W{rf-rf,hf)]AV, 



(5) 



NPe 

DtEyirf) = i E [HArf)D,W{rf - rf, hf)- 



H,{rf)D,W{rf - rf, hf)]AVj 



NPe 

DtE.irf) = 1 E [Hy{rf)D^W{rf - rf, hf)- 
H,{rf)DyW{rf - rf , hf)]AV, 

where rf and rf denote the positions of the z-th and j-th particle for the elec- 
tric and magnetic field component, respectively; NPe {NPh) is the number 
of neighbouring particles of the i-th fixed electric (magnetic) particle. The 
first derivative operator is indicated as Dx,Dy,Dz, x, y,z as the spatial co- 
ordinates. A central finite difference algorithm is used to approximate time 
derivatives in (jl]) and (jS]). A leapfrog scheme requiring only explicit updates, 
has been used for advancing the solution in time: the electric field is com- 
puted at whole time step and the magnetic field at the half time step [l|. 
By indicating with At the time step and with n = nAt the generic step, a 
complete temporal step of ([1]) evolves from n-1/2 to n+1, by involving both 
electric and magnetic field updates. The marching-on in time is actually 
from the step n-1/2 to the step n+1/2 when updating magnetic field while 
is from step n to the step n+1 when updating electric field: 



NPh 

H-^"\rf) = H2-'/\rf) -f^ [mrf)DyW{rf - rf, hf) 






(6) 



E'^{rf)D,W{rf - rf, hf)]AVj 



NPh 

Hr"\rf) = Hr''\rf) - f ^ [E:^{rf)DM{rf - rf , hf) 

3 = 1 



E^{rf)D,W{rf - rf, hf)]AVj 



NPh 

Hr'^\rf) = Hr"\rf) - f E \E-{rf)DM{rf - rf, hf) 

i=i 
E^{rf)DyWirf - rf, hf)]AVj 



NPe 

Er\rf) = E-{rf) + f E [Hr'^\rf)DyW{rf - rf, hf)- 



NPe 

E;+\rf) = E;irf) + f E [Hr'^\rf)D^Wirf - rf, hf)- 
H:^^'\rf)D^W{rf - rf , /if )] AV,- 



NPe 

E-^\rf) = E-{rf) + f ^ [H^^'^'{rf)D,W{rf - rf, hf)- 

m^^'\rf)DyW{rf - rf, hf)]AV, 

As well-known, the explicit time integration scheme is subjected to the 
CFL condition, that for the EM problem is expressed as: 

At< ^ 

cvrf 

where c is the propagation velocity of the physical phenomenon in the medium, 
d is the geometric dimension of the problem, and Ar is the minimum points 
spacing in the problem domain. In SPEM, the previous constraint depends 
on the smoothing length value; moreover, this last is strictly related to the 
spatial particles distances. Therefore it requires the time step to be propor- 
tional to the smallest spatial point resolution: At < mini{hi/c). 

The previous condition guarantees the numerical stability, but it can 
makes the SPEM impractical to use, especially when the distance between 
the two closest nodes becomes very small. In this section, some efforts are 
made to obtain a better maximum allowable time step. Namely, referring to 
the first case study reported in section 4, a propagation of a 2D TM^ wave 
is considered. In this case, the discrete formulation of relations ([6]) and (jTj) 
are compactly re-written, as follows: 

gn+i = E!? + At[TH"+^/2 _ uH^+i/2] (8) 



Jjn+l/2 _ jjn-1/2 _ ^^y^n ^g) 



H«+i/2 = u^-y^ + AtZE^ (10) 



where E^, H^:, H^ are the vectors of the field components and T, U, V, Z as 
the matrices collecting the spatial field derivatives computed with SPEM 
formulation and the physical parameters of the problem domain. By substi- 
tuting dl]) and ([TO]) in ®, i.e.: 



E;^+^ = [I + At2(TZ + UV)]E^ + At[TH;;~i/2 - UH^-^/2] (11) 

the evolution in time is described by the relation X""*"^ = SX"", where: 

I + At2(TZ + UV) -AtU AtT 
-AtV I 

AtZ I 

/ E^i \ / e: 

xn+i ^ H^+1/2 , , x'^ = nr^/^ 

In order to X"^^ be bounded, as n — > oo, the time increment must satisfy 
the following relation [27]: 



At < 



\/ >^MAxiO) 



:i2) 



where Amax(S) is the maximum eigenvalue of the matrix S. The flT2l) gives 
an upper limit to the increment of the time step getting away from the 
minimum particles spacing. Furthermore, the restriction on At limits the 
physical applications yet, so that in the following an unconditionally stable 
time advance equations has been introduced. The original SPEM method 
has been modified by preserving the leapfrog evolution in time of the explicit 
time integration scheme. 



3. The leapfrog ADI-FDTD method 

In this section the leapfrog ADI-FDTD method [Sj is generally described 
and re-formulated for generating the implicit leapfrog formulation of SPEM. 
In the ADI-FDTD method, the temporal stepping from n-1 to n step is 
obtained by splitting the integer time step into two half-time steps. By 
indicating with E and H the fields components with respect to a generic 
coordinate in the x,y,z geometric system, the finite difference time domain 
scheme involves the following general relations, in which the actual field value 
is formally expressed as a generic function of the fields at previous time steps: 



^-1/2 _ f(^E^-\ ir-i/2, ir-^) (13) 



ir-^/^ = g{ir~\Br~^/^,Er^^) (i4) 



BT" = /(£"-^/^ iT, ip-y^) (15) 



IP = g{ir~^/\F',E^-^/^) (16) 



For the marching-on in time of the electric field components, equation 
f lTB]) is substituted in f lT^ and fITB]) is substituted in flT^ . so obtaining: 

BT-^/^ = (t){Br-\H'-^) (17) 



In flT7|) and flTSj) first and second order spatial field derivatives hold and 
tridiagonal linear systems have to be solved to generate £""^" and E"^. An 
ADI-FDTD time step is so described by means of relations flTTl) and flT^ . 
and relation ()18j) and fITB]) . respectively. By considering the temporal step 

9 



for the fields evolution from step n — 1/2 to n and obtaining H" — 1/2 from 
(ITB]) . the following relation holds: 

ir-i/2 =g(i/^, £^, £"-i/2). (19) 

By substituting flT9l) into flTSl) the electric field component is generated 
by means of: 



E^ = tP{E"-^/^,H") (20) 

Let now consider equation (fT7|) for the marching-on in time of the electric 
field from n to n + 1/2: 



£"+1/2 = 0(£r,,^ ^) 

and by using fl20l) . the previous relation can be re- written as: 

r^'^' = 4>{E'-'/Mr) (21) 

The leapfrog equation describing the marching-on in time of the electric field 
component is so obtained. Let now consider relation ( ITSl) for the half time 
step from n-|-l/2 to n+1: 

rn+l _ ff pvi+l/2 IJri+l Tjn+l/2\ 

By substituting it into fITB]) . the following relation holds 

ii'^^ = g{W+^l\ Br+\ £"+V2) ^ ^(if'+l/2^ £n+l/2) (22) 

For the half time step from n to n+1/2 relation flT3l) is considered: 
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By substituting this last into flT4|) . the _f/"+i/2 component holds: 

ir+i/2 _ g(^^^ ^+1/2^ /(£"+'/', ir+1/2, iT)) (23) 

This last relation is subsituted in (|23|) so obtaining the following: 

ir+^=7(ir,£"+i/2). (24) 



Equations fl2T]) and (1251) give rise to LAF scheme for marching-on in time 
of the generic magnetic field component H. In (!2T|) and (l23l) only first and 
second spatial field derivatives are involved. The electric and magnetic fields 
are staggered in time and updated in one full time step, iteratively. 

3.1. The leapfrog ADI-FDTD method m SPEM 

The complete formulation of Maxwell's curl equations derived from eqs 
fl2T|) and flM|) is reported in the following by approximating the spatial first 
and second derivatives by means of SPH particle approximation: 

NPe 

eE:^''\rf) - (f )2i E E:^'/\rf)DlW{pf^hf)^V, = 

i=i 

NPh 

eE:-'/\rf) - (f )2i Y. E-.^"\rf)DlW{pf, hf)AV,+ (25) 

NPh 

At J2 [H^{rf)DyW{pf,hf)-H;{rf)DM{pf,hf)]AVj 



eE-^"\rf) - (f )^i E Er"\rf)DlW{pf^hf)AV, = 

i=i 

NPe 

eEy-"\rf) - (f )2i E Er''\rf)DlW{pf, hf)AV,+ (26) 

NPh 

At E [H-{rf)D,W{pf,hf) - H^{rf)D,W{pf,hf)]AV, 
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NPe 






NPe 

eE--"\rf) - (f )2i Y. Er"\rf)DlW{pf, hf)AV,+ (27) 

i=i 

NPh 

NPh 

eiJri(rf ) - (f )2i E ^r^rf )Z)^W^(gf , /^f ) AV, = 

NP„ 

^H^irf) - (f )2i E H-{rf)DlW{qf, hf)AV,+ (28) 

f E [i?r'/'(rf)/^.l^(gf,/^f) -i?r^/'(rf)D,l^(gf,/^f)]AK, 

e^;(rf ) - (f )^^ E ^,nrf )D,W(gf , hf)AV,+ (29) 

i=i 

e^r^(rf ) - (f )^^ E ^r^(rf)/^^w^(gf,/^f)A\/,= 

AfPa 

ei^."(rf ) - {f?l E ^r(rf )I^^l^(gf , /^f )A\/,+ (30) 

A^Pb 

f E [E-.^"\rf)D,W{qfX)-ET''\rf)DM{qf,hf)]AV, 

where p|' = rf — r|', p|^ = rf — r^ , if = rf — rj, q^ = rf — r^ . 

The equations (25)- (30) can be expressed in a compact matrix notation as 

follows: 

where: 
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\ 











-^yy 











K. 


^xx 











^yy 











c. 



A = A„„ Lb 



c 






^xy 


B.. 


Bya; 





-B,. 


^zx 


Bzy 








— F 


F., 


F.. 





-F,. 


— F 

-■- zx 


F,„ 






and the non zero entries are: 



I3 = y,z,x I b,, (t, j) = k, (^, j) = AtD.W {pf, hf) AV, ^ J^^y^'^^^^^j, 



The matrices employed in A,B,C,F are sparse block matrices and the 
sparsity depends on the influence domain of each particle. Each row of 
the matrices refers to a fixed particle, center of a kernel function, and the 
non zero elements correspond to the particles surrounding the fixed one. 
The matrix A^x {-^yy, A^^) is generated from the smoothing kernel functions 
regarding the E'-particles , i.e. each kernel is centered on rf and has Tj" as 
neighbourings. In a similar fashion the matrix Cxx {Cyy,Czz) is related to 
the smoothing kernel functions centered on r^ and has Vj^ as neighbourings. 
In all these matrices the second derivatives are involved. The matrices Bj^./ 
{k = x,y,z;l = x, y,z,l ^ k) regard smoothing kernel centered on rf but 
involving r^ as neighbourings. Analougously, F^; {k = x,y,z;l = x, y,z,l j^ 
k) regards smoothing kernel fixed on rf but involving r^ as neighbourings. In 
all these last matrices the first derivatives are involved. When kernel vanish, 
the limit configuration of the explicit formulation in time is obtained. In the 
following, some matrices configurations are reported. Equally distributed 
particles are considered generating a structured block matrices (Figs. 1,2). In 
this case the matrices are banded and the bandwidth is with the width of 
the local support of the kernel function. The smoothing length h=aAr must 
be opportunely chosen to achieve satisfactory approximation field values, 
giving rise to a sparse and solvable system. The approach in which the 
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support size remain fixed for increasingly denser set of data is very inefficient 
witli associate system matrices increasingly denser. In the simulations, the 
smoothing length is scaled so that peaked basis functions are with densely 
spaced particles and flat basis functions are with coarsely spaced particles. 
In Fig. 3 the error in the computation is reported by considering 100 equally 
distributed points over the whole problem domain. For the electric fleld in the 
TMz case, the error is measured in L2 in comparison with the flnite difference 
time domain solver. The best behaviour is with q;=0.7075. In tab.l results 
are reported by considering the optimal smoothing length found, to give an 
accurate solution for this problem. The sparsity of the matrices increases by 
increasing the problem size. 



N 


Sparsity 


100 


81.86 % 


200 


91.412 % 


625 


96.45 % 


900 


96.45 % 


2500 


97.48% 


3600 


99.33% 



Table 1. Results with regular particle distribution by increasing the parti- 
cles density N over the problem domain in the computation of the matrix A^x 

In Figs. 4, 5 the block matrix A^x is depicted for irregular particle distri- 
butions, by varying the number of particle in the problem domain. Fig.4 is 
with 81 particles and a sparsity of about 80 % is generated. Fig.5 refers to 
625 particles with sparse matrix of about 96 % of zero elements. 

The kernel derivatives in computing the derivatives of the fleld compo- 
nents can be performed in a pre-processing stage, and all the matrices em- 
ployed in the computation can be assembled out of the temporal loop. More- 
over, by considering the smoothing length independent from the evolution in 
time, the structure of the matrix not depends on the temporal step yet. 

4. Numerical validation 

In this section, some different case studies have been performed and re- 
sults are discussed to assess the proposed meshless implicit leapfrog numerical 
method. At flrst, a transverse electric (TM^)) mode in air is flrstly taken 
into account. In this case, equations are reduced as follows: 
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Figure 1: Profile of the matrix A for a regular particle distribution 
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Figure 2: Zoom of the skeleton of the matrix Axx in A of figure 1 
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Figure 3: Error behaviour by varying h for 100 equally distributed particles over the 
problem domain. 



DfEz = -{D^Hy — DyHx), DtHy = -{DxEz), DtH^ = -{DyEz 



(32) 



A regular particles distribution is built up, such as in a FDTD spatial 
grid. For problem with regular distribution, nodes are placed as in the 
point-matching procedure; for the general case, a Voronoi decomposition 
[2 8(1 or others technique can be adopted. In particular, in the paper when an 
irregular particles distribution is set, in order to avoid particle inconsistency 
which can lead to a loss of accuracy the consistency restoring adopted in 
pj, |2|, is used. The FDTD simulation is considered for comparison. The per- 
fect matching layer (PML) technique is used in simulating the open spatial 



domain [ly, |28| . In this way, also a good treatment of the surface integral 
truncation problem, as reported in section 2, is carried out for LAF-SPEM 
computation. The actual time step, satisfies the CFL condition, so as in the 
FDTD approach, i.e.: 



At = AcFL 



Ax 
2c 



(33) 
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Figure 4: Profile of the matrix A^^^: for an irregular 9x9 particle distribution 
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Figure 5: Profile of the matrix Kxx for an irregular 25x25 particle distribution 
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where Ax = Ay = 1cm is the spatial step. In figure 6 FDTD and LAF- 
SPEM spatial profiles along y direction, ten cells distant from the middle in 
the X direction in a 2D domain, 100x100 number of cells are reported; a time 

(20-t)^ 

varying shaped source, e 72 , placed at the middle point of the domain is 
considered in simulating a TM^ field in air. The comparison shows a satis- 
factory agreement. 



titre-step n. 70 




Figure 6: Comparison among FDTD and LAF-SPEM results in a 2-D domain 100x100 
cells and particles: the time dependent source is placed at the middle point and a TMz 
field is simulated in air. E^ field (V/m) spatial profiles at time step 70, along y direction 
ten cells distant from the middle in the x direction 



As a further validation, an axial symmetric cylindrical domain in air is 
considered, with Tq = 0.2m (eo,/io as constitutive parameters), with the 
following boundary and initial conditions: 



a/x^ + y^, To = 0.2m, < r < rg 



(34) 



18 



E,{x,y,0) 



'o 



n SE4x,y,0) _ p, 
^' 5t ~ ^ 



(35) 



The domain is discretized with 10x10 irregular distributed particles ob- 
tained by randomly positioning the particles near the original regular position 
in the cylindrical domain. The actual time step is set equal to AAtcFi- As 
shown in figure 7 the simulation remains stable also for a final time cor- 
responding to 175 time steps, and a satisfactory approximation has been 
achieved. 



time step n. 175; final time=4.6688e-008 [s] 




10 15 20 25 30 35 40 45 50 



Figure 7: Cylindrical domain with tq — 0.2m bounded by a perfect electric conductor: 
comparison between LAF-SPEM results and the exact solution at time step 175. The 
LAF-SPEM time step is 4 times that satisfying the CFL condition. 

In order to better assess the validity of the proposed approach, a more 
realistic example has been carried out. A sectorial (2-D) Perfectly Elec- 
tric Conducting (PEC) horn antenna excited by a sinusoidal voltage in a 
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transverse magnetic (TM^) computational domain has been simulated and 
compared with FDTD results [16] [lO|. The computational domain is trun- 
cated with a PML absorbing boundary condition. The horn is modelled by 
a perfect electric conductor (PEC) material. The spatial cell size is equal to 
0.0025 m, the time step is equal to 4.23 ps, the excitation frequency is equal 
to 9.84 GHz. 



LAF-SPEM at time:8.382e-010 [s] 
L relative error=3.71 % 
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Figure 8: Sectorial (2-D) PEC horn antenna excited by a sinusoidal voltage in a 2-D 
domain with 100x100 regular particles distribution. Colour map of the Ez field (V/m) at 
time step f 98 obtained with LAF-SPEM - relation ([33|) holds. The reported relative error 
is related to FDTD results. 



In figure 8, LAF-SPEM E^ field results, in a 2-D domain with 100x100 
regular particles distribution, are reported, at time step 198. The excitation 
section is also reported in the figure. The actual time step. At , satisfies the 
CFL condition (l33l) . In comparison with FDTD simulation, a good agreement 
has been obtained. For a better comparison with FDTD results, in figure 9 
FDTD and LAF-SPEM Ey field spatial profiles along the direction containing 
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the longitudinal axis of the horn antenna at the same time step of figure 6, 
are reported. 
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Figure 9: Sectorial (2-D) PEC horn antenna excited by a sinusoidal voltage in a 2-D domain 
with 100x100 regular particles distribution. Comparison between FDTD and LAF-SPEM 
E field (V/m) results at time step 198, along the direction containing the longitudinal axis 
of the antenna - relation psp holds. 

In order to show the capability of the proposed approach to run with 
a time step greater than that obtained by satisfying the CFL condition, a 
simulation with a time step equal to 4 At CFL is carried out. Also in this 
case, an irregular particles distribution has been also adopted by randomly 
positioning the particles near the original regular position in the domain. 
As shown in figure 10, the simulation remains stable also for a final time 
corresponding to 200 time steps, and a satisfactory approximation has been 
achieved. 

5. Conclusions 

In this paper the SPH meshless method has been composed with a leapfrog 
ADl-FDTD method for time evolution of electromagnetic time-domain tran- 
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fral timE:3.3782e-009 [s] 




Figure 10: Sectorial (2-D) PEC horn antenna excited by a sinusoidal voltage in a 2-D 
domain with 100x100 irregular particles distribution. Comparison between E^ field (V/m) 
results obtained with FDTD and LAF-SPEM at the same final time, along the direction 
containing the longitudinal axis of the antenna. The LAF-SPEM time step is 4 times that 
satisfying the CFL condition. 



sient propagation problems. The new schemeimproves the original SPEM 
method, already developed by the authors. The method is very promising 
because avoids grid generation in space and results unconditionally stable 
in time. In this way, realistic simulations can be performed with reasonable 
computational efforts. The computational tool is assessed by using different 
numerical simulations also employing irregular particles distribution. Perfect 
matching layer technique is used in simulating open spatial problems and a 
consistency restoring approach is introduced. The proposed methods may 
have good perspectives of extensive applications. 
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